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Grid  Discretization  Based  Method  for  Anisotropic  Shortest  Path 
Problem  over  Continuous  Regions 

Zhanfeng  Jia  and  Pravin  Varaiya 


Abstract — The  paper  presents  a  method  to  find  the  shortest 
path  between  two  points  over  a  continuous  region.  The  length 
of  a  path  is  the  integral  of  the  cost  along  the  path.  The 
cost  can  be  anisotropic,  meaning  that  It  depends  on  both  the 
position  on  the  path  and  its  direction.  The  method  uses  a 
simple  rectangular  grid  to  discretize  the  region,  independent 
of  die  cost  Because  the  cost  is  anisotropic,  rectlinear  paths 
connecting  adjacent  grid  points  may  not  approximate  the 
optimal  path.  To  overcome  this  ‘iimit-on-directlon’’  problem, 
the  method  searches  over  shifted  versions  of  the  rectlinear 
paths.  A  Bellman>Ford  style  algorithm  finds  the  best  shifted 
path.  Theoretical  analysis  and  numerical  experiments  ensure 
efficiency  of  the  algorithm. 

1.  Introduction 

The  problem  is  to  find  the  shortest  path  over  region  71 C 
For  a  starting  point  €  7^,  the  ‘length’  V(x,p)  of  a 
continuous  path  p  from  Xs  to  a  destination  x  is  the  integral 
of  the  cost  function  c(x,  u)  along  the  path  p, 

V{x,p)^  r  c{p{t),p{t))dt.  (1) 

Jo 

The  path  p  :  [0,T]  T^^liiEst  b¥  contmubuTly^Wice 
differentiable;  p(0)  =  x,, p(T)  =  x-^and  such  that  |lj5||  =  1, 
and  IIpII  is  bounded.  We  wanra'mbfhod  fflSf  finds Ih?  best 
path, 

V’(x)  =  min  ^(xjp).  (2) 

peP(x,,x) 

V{x)  is  the  value  function  at  x.  The  cost  c  :  71  x 
(0,  oo)  depends  on  both  position  p{t)  and  velocity  p. 
However,  as  ||p||  =  1,  the  cost  c  depends  only  on  position 
and  direction.  Problem  (2)  is  called  the  Anisotropic  Shortest 
Path  (ASP)  problem.  The  special  case,  when  c(p,p)  is 
independent  of  p,  is  the  Isotropic  Shortest  Path  (ISP) 
problem. 

The  ASP  problem  arises  naturally  in  path  planning.  Our 
motivation  is  to  find  paths  for  Unmanned  Aerial  Vehicles 
(UAVs)  that  can  reach  a  designated  target  while  minimizing 
the  risk  of  detection  from  enemy  radar.  The  radar  signature 
of  a  UAV  is  a  function  of  its  distance  from  the  radar  and 
its  flight  angle,  so  the  accumulated  risk  faced  by  the  UAV 
along  a  path  can  be  expressed  as  in  (1). 

One  approach  to  ASP  problem  is  to  select  a  uniform 
grid  of  points  spaced  e  apart,  with  a  link  between  any  two 
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nearby  grid  points.  The  cost  of  a  link  connecting  points  x,  y 
is  taken  to  be 

c(x,irf^)x||x-j,||,  (3) 

II  a;  -  y  II 

and  one  can  use  an  algorithm  like  Dijkstra’s  to  find  the 
shortest  path  across  fhe  graph.  Kim  and  Hespanhajl] 
show  that  the  shortest  path  along  the  graph  convei^es  to 
a ‘shortest  continuous  path  as  the  graph  is  refined.  For  a 
two-dimensional  region,  the  graph  constructed  above  has 
0(e“')  vertices  and  0(6“^)  links,  which  is  inefficient. 

Kim  and  Hespanha  [1]  propose  an  efficient  adaptive 
sampling  scheme  that  takes  into  account  the  characteristics 
of  the  cost  function  over  the  region.  Using  a  uniform  grid 
is  inefficient  because  of  the  “limit-on-direction”  problem, 
discussed  in  the  next  section. 

The  proposed  method  finds  approximately  shortest  paths 
on  a  graph  with  a  uniform  grid.  The  crucial  difference 
from  the  graph  described  above  is  that  there  is  a  link  only 
between  adjacent  vertices.  Consequently,  the  graph  is  easy 
to  construct,  and  unlike  in  (3)  there  is  no  need  to  evaluate 
the  distance  between  pairs  of  point  in  order  to  select  the 
links  in  a  shortest  path.  The  algorithm,  called  the  Shiftable 
Segment-by-SegmerU  Path  (SSSP)  algorithm,  overcomes  the 
limit-on-direction  problem  by  allowing  the  paths  to  shift 
away  from  the  grid  points.  The  shifts  are  local,  so  the 
shifted  path  is  a  close  variant  to  the  discrete  path,  and  better 
approximates  the  desired  shortest  path. 

Section  II  provides  a  brief  review.  Section  III  introduces 
the  segment-shifting  idea  and  discusses  how  it  helps  to  find 
the  shortest  continuous  path..  Section  IV  describes  the  SSSP 
algorithm  and  analyzes  its  complexity.  Section  V  presents 
numerical  experiments.  Section  VI  concludes  the  paper. 

II.  Isotropic  and  Anisotropic  Shortest  Path 
Problems:  Related  Work 

An  efficient  Dijkstra-like  algorithm  for  the  Isotropic 
Shortest  Path  (ISP)  problem  was  first  proposed  by  Tsitsildis 
[4]  in  1995,  and  later  independently  proposed  by  Adal- 
steinsson  and  Sethian  [2]  and  called  the  Fast  Marching 
algorithm.  Using  a  uniform  grid,  the  algorithm  mimics  a 
wave  front  propagation.  The  idea  comes  from  geometric 
optics  where  a  wave  front  propagates  at  different  speeds 
depending  on  the  medium  which,  in  this  case,  is  inversely 
proportional  to  the  cost  function  at  each  point  x.The  cost 
does  not  depend  on  the  direction  p(i),  so  Problem  (2) 
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Fig.  1.  Value  function  update  in  Ae  Fast  Marching  meAod:  V{x)  — 
V{z)  +  c(x)  •  Ien(x2),  where  V{z)  is  intcrylated  by  Ae  two  adjacent 
pomts  of  X,  and  len(x2:)  = 


Fig.  2.  Failure  of  Ae  FM  method  on  an  anisotroinc  case  wiA  elliptic 
cost  profile.  While  Ae  shortest  paA  comes  from  between  point  j4'  and 
the  result  of  FM  is  between  pomt  B  and  D. 


simplifies  to 

V{x)  =  min  f  c{p{t))dt.  (4) 

p€P(Xa  ,*)  Jq 

V{x)  can  be  viewed  as  the  time  the  wave  first  arrives  at 
X.  TTie  Fast  Marching  algorithm  strategy  uses  the  value 
function  to  order  the  selection  of  grid  points,  similar  to 
Dijkstra’s  algorithm,  as  the  current  smallest  wave  front  The 
underlying  theory  is  the  causality  relationship^  which  states 
that  the  arrival  time  (value  function)  V{x)  at  x  depends 
only  on  the  neighbors  with  strictly  smaller  F-values.  The 
proof  of  this  statement  can  be  found  in  [4]* 

Using  a  uniform  grid  and  the  causality  relationship,  the 
Fast  Marching  method  updates  the  value  function  at  a  grid 
point  x  according  only  to  the  values  of  its  adjacent  grid 
points.  Specifically,  in  the  two  dimensional  case  the  update 
process  is  given  by 


V(x)  =  min  min 

'  a€A0<3<t 


4-  +  haici)  -f  (1  —  s)V{x  -f  ha2e2) 


,  (5) 


where  ei  and  62  are  the  unit  vectors  in  a  =  (ai,a2) 
is  an  element  of  A  =  {-1,1}^,  and  h  represents  the 
discretization  step  (see  Fig.  1).  Equation  (5)  resembles  a 
discrete  approximation  of  the  Hamilton-Jacobi  (HJ)  equa¬ 
tion. 

The  Fast  Marching  algorithm,  however,  fails  in  the 
anisotropic  case,  one  reason  being  the  failure  of  the  causal¬ 
ity  relationship  on  the  grid,  as  illustrated  by  Sethian  and 
Vladimirsky  [3].  In  their  example,  the  cost  function  has 
an  elliptic  profile  and  is  uniform  across  the  region.  The 
constant-value  contours  are  concentric  ellipses  with  center 
Xa,  and  the  shortest  path  to  any  point  x  is  the  straight  line 
connecting  Xs  and  x. 

The  Fast  Marching  algorithm  cannot  converge  to  these 
shortest  paths.  Suppose  we  want  to  update  the  value  func¬ 
tion  at  point  C  in  Fig.  2.  The  thick  elliptic  curve  is  a 
constant- value  contour,  so  V (A)  is  larger  tfian  V (C),  and 
both  V(B)  and  V{D)  are  less  than  V{C).  Given  the 
causality  relationship,  the  Fast  Marching  algorithm  uses 
points  B  and  D  to  update  V{C)  and  omits  the  combination 
of  points  A  and  B.  Thus,  the  path  suggested  by  the 
Fast  Marching  algorithm  can  not  approximate  the  optimal 
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Hg.  3.  Limit-cm-direction  problem  of  Ae  grid  sampling.  For  a  grid 
discretization  of  step  h,  and  a  neighborhood  distance  S  =  y/2h,  there 
are  eight  neighbor  points  of  x,  implying  only  eight  possible  directions. 
The  dashed  circle  illustrates  the  neighboAood  area  of  point  x. 


continuous  path.  Modifications  of  the  Fast  Marching  method 
have  been  proposed  to  address  this  problem.  For  example, 
in  [3],  the  update  procedure  considers  not  only  the  four 
adjacent  points,  but  non-adjacent  points  such  as  A'  and  B\ 
at  the  cost  of  larger  time  complexity. 

Another  difficulty  introduced  by  anisotropy  comes  from 
the  restriction  that  ||p||  is  bounded.  The  restriction  is 
reasonable  in  the  UAV  scenario,  in  which  the  UAV’s  ability 
to  change  directions  is  limited.  But  the  bound  leads  to  a 
failure  of  the  Principle  of  Optimality^  because  a  partial 
section  of  a  shortest  path  may  not  be  itself  be  a  shortest 
path.  Therefore,  the  value  function  V(x),  as  a  function  of 
position  alone  does  not  propagate  properly.  One  can  extend 
the  value  function,  with  bofii  position  and  direction  as  state, 
but  this  Fast  Marching  much  more  complicated. 

Alternative  approaches  to  the  ASP  problem  focus  on 
paths  instead  of  value  functions.  The  idea  is,  first,  to  reduce 
the  continuous  region  to  a  discrete  graph,  and  then  to 
compute  the  shortest  paths  across  the  graph.  Path-oriented 
methods  fit  the  ASP  Problem  (2)  well  because  paths  cany 
direction  information.  Moreover,  path-oriented  methods  re¬ 
move  the  additional  procedure  needed  to  reconstruct  paths 
from  the  value  functions. 

The  basic  step  associated  with  path-oriented  methods  is 
sampling,  which  reduces  the  region  into  a  discrete  graph. 
The  simplest  is  (uniform)  grid  sampling.  Unfortunately,  grid 
sampling  is  not  efficient  in  the  sense  that  the  number  of 
vertices  and  edges  grows  very  fast  if  we  want  the  shortest 
discrete  path  to  approximate  both  the  position  and  direction 
of  the  shortest  continuous  path.  This  is  because  of  the 
“limit-on-direction”  effect. 

Remember  that  the  discrete  paths  are  composed  of 
straight  segments  between  sample  points.  Consider  a  two- 
dimensional  grid  and  a  small  neighborhood  distance  6  = 


2^ 


V2h  as  in  Fig.  3.  The  dashed  circle  centered  at  x  shows 
the  neighborhood  of  x.  Denote  the  eight  grid  points  in  the 
neighborhood  by  Suppose  only  these 

eight  points  in  are  available  to  form  a  segment  with 
X,  Then  every  possible  segment  in  the  grid  graph  has  to  lie 
either  along  grid  lines  (as  from  2/2*  2/4*  Vs*  or  yr  to  x),  or  on 
diagonal  lines  of  grids,  (as  from  yi,  2/3,  ye,  or  ys  to  x),  thus 
on  one  of  the  eight  directions.  The  discrete  paths  composed 
of  these  segments  can  never  approximate  a  continuous  path 
along  (say)  a  30  degree  line.  Increasing  the  fineness  of 
discretization  or  the  neighborhood  distance  both  introduce 
noore  directions,  but  cannot  eliminate  the  limit-on-direction 
effect  completely.  In  fact,  even  the  random  sampling  works 
more  efficiently  than  grids. 

Kim  and  Hespanha  in  [1]  proposed  a  non-uniform  sam¬ 
pling  scheme  called  Honeycomb,  which  allows  sparse  sam¬ 
pling  at  areas  where  the  gradients  of  the  cost  function 
over  position  and  direction  are  both  small.  The  Honeycomb 
scheme  first  generates  cells  whose  diameters  are  inversely 
proportional  to  a  linear  function  of  the  two  gradients. 
The  sampling  points  are  then  put  on  the  cell  border.  The 
efficiency  of  the  Honeycomb  scheme  comes  from  its  careful 
dependency  on  the  cost  function. 


III.  Shifted  Paths 

We  modify  the  discrete  path  based  on  grid  discretization 
and  allow  the  end  points  of  the  segments  to  shift  away 
from  the  grid  points.  The  shifts  are  local,  so  the  shifted 
path  is  a  close  variant  to  the  original  discrete  path.  With  a 
shift  the  segments  can  take  virtually  any  directions,  thereby 
overcoming  the  limit-on-direction  problem  associated  with 
grid  sampling. 

Recall  the  formulation  of  the  path-oriented  method  in 
[1].  The  compact  region  is  sampled  by  a  finite  set  of 
points  X  which  includes  the  starting  point  x^.  The  discrete 
path  p(t)  over  X  is  defined  to  be  a  piecewise  linear  curve 
connecting  a  sequence  of  points  in  set  X  with  respect  to 
a  given  neighborhood  distance  <5.  Specifically,  there  is  a 
sequence  of  points 

{xo  =  Xs,Xi,X2,...,X,n  =  x}  C  X  (6) 

that  satisfy  the  neighborhood  condition,  as  in  (7)  below, 
llx,^Xk-ill<S,Vk€{l,2,,..,m}.  (7) 

The  discrete  path  p(t)  is  defined  as 


p(t)  =  Xfc-i  -h 


Xjt  -  Xjt-l 


Ijxife  “Xjfc^ll 

Vt  €  [tk^uhlke  {l,2,...,Tn},  (8) 


where  tk  is  defined  recursively  by 

to  =  0,4  =  4-1  +  |jxjt-Xfc-i||,  Vk  €  (9) 

and  tm  =  T.  The  fcth  segment  of  p(t),  denoted  by 
or  TTfc  if  there  is  no  confusion,  is  the  part  of  path  p(t) 
with  t  €  (4,4+i]>  k  e  {0, 1}.  Points  Xk  and 


Xfc+i  are  called  the  start  point  and  end  point  of  segment 
iTk  respectively.  Note  that  ||p(f)  =  1||  is  assumed  in  the 
above  definition,  time  is  thus  not  important  in  most  cases. 
It  is  then  convinient  to  represent  the  discrete  path  p  as 
P  =  {xo  ==X5,xi,...,Xm  =  a:}. 

Henceforth  we  use  grid  sampling  with  grid  size  h  to 
discretize  the  region.  The  discrete  set  X  contains  all  points 
X  €  71  of  the  form  x  =  (xs,i  -f  4^,  +  inh)^ 

where  x,  is  the  starting  point  and  4,— ,4  arc  integers. 
The  shifted  path  can  now  be  defined,  starting  with  a  discrete 
path  p={xo  —  X5,xi,  ...,x,n},  in  which  the  x^’s  are  grid 
points  satisfying  (7).  A  shifted  path  p  with  respect  to  p  is 
defined  by  the  shift  array  5  =  (si,  ...,5m)  such  that 

P5  =  {^0,Xi,...,Xm}  =  {z:5,Xi-{-5i,...,Xm-|-5m},  (10) 

5jfe  €  ft  are  local  shifts,  and  Xk  €  ft(xfc)  :=  {xj^+s  :  s  6  ft} 
is  within  the  local  shift  area  of  xjt.  One  can  apply  definition 
(8)  to  obtain  the  shifted  path  ps{t)  as 

Vf  e  (11) 

where  4  is  defined  recursively  by 
*0  =  0,ffc  =ffc_i+||ifc-Xfc_i||,  Vfc  e  {l,2,...,7n}.  (12) 

The  local  shift  area  ft  is  the  /i-ball, 

ft  :={5:5€MM|5||<ft}.  (13) 

The  size  of  ft  matches  the  grid  size  h,  assuring  the  locality 
of  the  shifts.  Moreover,  provided  that  the  neighborhood 
distance  6  is  reasonably  large,  a  segment  of  the  shifted  paths 
is  able  to  take  any  directions  as  the  end  point  of  the  segment 
shifts  within  the  area  of  ft,  thus  overcoming  the  limit-on- 
direction  problem.  This  property  is  listed  as  the  following 
proposition. 

Proposition:  Suppose  the  neighborhood  distance  S  > 
y/nh  where  n  is  the  dimension  of  the  region  and  h  is  the 
discretization  step.  Consider  a  segment  tt  with  the  start  point 
X  ^  X  and  an  arbitrary  direction  v  {v  £  R”,  ||v||  =  1  spans 
the  surface  of  the  unit  ball).  There  exists  a  neighbor  point 
2/  G  34  and  a  shift  s  e  ft  such  that  the  end  point  of  segment 
TT  becomes  p  =  p  +  s,  and  the  direction  of  tt  is  same  as  the 
direction  v.  Specifically,  there  exists  a  scalar  ot  such  that 

y  X  =  av.  (14) 

We  omit  the  technical  proof.  The  proposition  is  intuitively 
obvious  if  y  spans  a  set  that  surrounds  x  completely,  as 
does  the  surface  of  the  unit  ball  centered  at  x.  Fortunately, 
given  5  >  y/nh,  34  contains  all  vertices  of  the  cube  of 
size  h  cornered  at  x.  With  ft  defined  in  (13),  the  set  {y  = 
y-l-s:ye34,s€ft}  contains  the  surface  of  this  /i-cube, 
therefore  surrounds  x  completely,  and  thereby  eliminates 
the  limit-on-direction  problem  of  the  grid  discretization. 
Though  the  definition  of  ft  in  (13)  is  sufficient  for  this 
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purpose,  the  set  may  be  too  large.  We  can  take  as  the  low¬ 
dimensional  set  that  restricts  the  shift  5  to  be  perpendicular 
to  the  path, 

n  :=  {s  :  s  €  ||5||  <  h,  {s,p)  =  0},  (15) 

and  still  overcome  the  limit-on-direction  problem.  Taking  Q 
as  in  (15)  reduces  the  number  of  parameters  associated  with 
a  shift  path  so  that  the  optimization  procedure  described 
below  remains  simple. 

The  concept  of  shifted  paths  gives  a  new  meaning  to 
the  value  function  of  a  discrete  path  in  the  context  of 
grid  discretization.  Given  a  discrete  path  p  ==  {xo  = 
=  a^},  the  tilde  value  function  V{x,p)  is 

defined  by 

V{x,p)^  V{x^Ps)-  (16) 

The  minimization  is  over  all  local  shifts  of  the  discrete  path 
p,  so  V (x,  p)  is  the  lowest  value  one  can  expect  by  choosing 
p.  Therefore,  even  when  the  grid  is  coarse,  the  minimizing 
shifted  path  p  has  a  good  chance  to  closely  approximate 
a  shortest  continuous  path.  As  we  refine  the  grid,  the 
optimal  V{Xfp)  will  converge  to  the  optimal  solution  of 
Problem  (2). 

Theorem:  For  any  A,  e  >  0,  there  exists  a  finite  set  X  of 
iVc  grid  sampled  points  of  7^,  due  to  the  discretization  step 
ft,  such  that  for  any  given  starting  point  Xs  ^  X  and  any 

X  ^ 

min  K(x,p)-€  <  ^(x)  <  min  V"(x,p),  (17) 
pePx  (®.  ,x)  p^Px  ,«) 

where  Px{xs,x)  contains  all  the  discrete  paths  from  Xs  to 
X  over  Xj  A  is  the  bound  on  the  second  derivative  of  the 
continuous  path  p,  and  V  (x)  is  the  value  function  of  the 
shortest  path  of  the  Anisotropic  Shortest  Path  problem,  as 
defined  in  (2). 

The  theorem  is  proved  along  along  the  same  lines  as 
Theorem  1  in  [1].  The  latter  guarantees  the  existence  of  a 
finite  set  X  of  size  such  that  one  can  construct  a  discrete 
path  po  over  X  satisfying 

V{x,po)-V{x)<€.  (18) 

Notice  that  the  definition  of  the  tilde  value  function  in  (16) 
implies  K(x,po)  <  F(x,po),  so  the  convergence  theorem 
holds  immediately.  However,  this  argument  doesn’t  reveal 
how  fast  the  num^r  of  points  increases,  or  equivalently, 
how  fast  ft  decreases,  as  e  decreases.  In  the  numerical 
experiments  the  shifted  paths  perform  much  better  than  the 
discrete  paths,  meaning  that  the  proposed  algorithm  can 
work  over  a  coarse  grid  with  a  small  number  of  points  in 
A',  and  yet  yield  paths  with  values  close  to  the  optimum. 

IV.  Algorithm 

The  underlying  scheme  of  the  proposed  Shiftable 
Segment-by-Segment  Pathing  (SSSP)  algorithm  is  the 
Bellman-Ford  architechture,  in  which  each  point  records 
the  most  promising  paths  and  their  tilde  value  functions 


Rg.  4.  Some  cxan^les  of  segments  (thick  lines  with  anow)  and  their 
baseline  directions  (below)  when  I  =  3.  is  the  shift  al<mg  the  grid  line 
that  is  peipendicular  to  the  baseline  direction.  The  dashed  lines  show  the 
shifted  segments. 

from  the  starting  point  and  tries  to  extend  them  to  the 
destination  through  its  neighbors.  The  paths  eventually 
reach  the  destination  and  the  shortest  path  (with  the  smallest 
tilde  value)  among  them  is  the  solution.  Note  that  the 
recorded  paths  are  discrete  paths.  To  find  the  tilde  values  of 
these  paths,  one  has  to  solve  an  optimization  problem  over 
the  shift  vector  Sk  for  all  segments  TTfc.  If  a  path  p  has  m 
segments,  it  is  an  optimization  problem  over  which  is 
obviously  complicated  and  time  consuming.  We  introduce  a 
heuristic  meth^  to  compute  the  tilde  values  approximately. 
From  now  on  we  restrict  the  region  to  a  two-dimensional 
plane.  According  to  (15),  the  shifts  s^’s  are  now  scalars. 

Consider  the  graph  formed  by  a  rectangular  grid  with 
only  vertical  and  horizontal  links  between  adjacent  points. 
Let  the  segment  tt  be  a  group  of  links  of  up  to  I  base  links 
and  0  or  1  perpendicular  link,  as  shown  in  Fig.  4.  When 
/  =  1,  it  corresponds  to  the  case  where  the  neighborhood 
distance  6  =  v^ft;  and  1  =  3  implies  S  =  y/l6h.  The 
direction  of  the  base  links  is  called  the  baseline  direction 
of  the  segment,  while  the  direction  the  perpendicular  link, 
if  there  is  one,  is  called  the  quadrantal  direction.  The  shift 
Sk  of  the  start  point  of  a  segment  tta*  is  along  the  grid  line 
that  is  perpendicular  to  the  baseline  direction  of  itk^  The 
assumption  here  is  that  the  baseline  direction  represents  the 
direction  of  the  path. 

The  main  loop  of  the  SSSP  algorithm  follows  the 
Bellman-Ford  architecture.  For  a  graph  of  n  =  \X\  points, 
the  architecture  does  up  to  n  -  1  rounds  over  all  links 
relaxing,  or  updating,  the  tilde  value  function  (the  distance 
to  the  starting  point  Xj),  as  well  as  the  shifted  path  that 
generates  it  This  information  forms  a  record  of  a  point 
X  e  A",  including  the  discrete  path  p* ,  the  shifts  S'*,  and  the 
corresponding  tilde  value  function  V*(x,p*).  In  practice, 
the  relaxation  operation  is  initiated  by  a  recently  updated 
point  Xjfc  when,  from  a  path-oriented  view,  its  recorded 
path  Pit  :=  Pxk  =  {xsyXi,.„,Xk^i)Xk}  tries  to  extend 
to  a  neighbor  point  xjt+i.  Denote  the  extended  path  by 
Pk^i  :=  {pife,Xit+i};  the  records  at  Xfc+i  will  be  relaxed, 
or  updated,  if  the  tilde  value  of  pk^i  is  smaller  than  the 
recorded  V*{xk^i,Px^^i)^  Initialized  to  infinity,  the  tilde 
value  function  therefore  decreases  with  each  successful 
relaxation  operation.  The  following  procedures  describe  the 
details  of  the  relaxation  operation  of  SSSP. 

1)  Approximate  the  tilde  value  function  V{zkiPk)  of  the 
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point  Zk  that  shifts  away  from  the  grid  point  Xk* 
Zk  =  Xife  +  is  the  shifted  point,  H  is 

the  shift,  and  is  the  unit  vector  representing  the 
direction  that  is  perpendicular  to  the  baseline  direction 
of  the  segment  Ttk  =  This  is  the  key 

computation  of  SSSP  and  we  will  discuss  it  in  detail 
soon. 

2)  Compute  the  tilde  value  function  K(xjt+i,pfc+i)  of 
the  extended  path  pit+i  as 

t^(a;fe+i,Pfc+i)  =  mm{V{zk,Pk) 

+  INfe+i  -  Zk\\  •  c(xfc+i, Xfc+i  -  Zk)),  (19) 

where  c(x,  v)  is  the  cost  function. 

3)  If  V{xk^i,pk^x)  <  V*{xk^upl^^J,  the  relaxation 
operation  succeeds  and  the  record  of  xjk+i  is  updated. 

The  Bellman-Ford  architecture  of  SSSP  allows  each  grid 
point  to  keep  one  record  of  the  most  promising  path.  How¬ 
ever,  because  of  the  failure  of  the  Principle  of  Optimality, 
a  sub-optimal  path  at  Xk  may  extend  to  be^e  best  path 
at  Xk^i>  Therefore,  it  is  efficient  to  keep  multiple  records 
of  the  most  promising  paths  at  each  points.  The  size  of  the 
records  doesn’t  need  to  be  large.  For  example,  when  the 
number  of  base  links  of  a  segment  is  selected  to  be  Z  =  1, 
there  are  at  most  eight  different  choices  of  neighbor  points. 
Therefore,  each  point  needs  to  keep  no  more  than  eight 
paths,  more  paths  make  no  difference.  Moreover,  keeping 
multiple  records  increases  the  time  complexity  of  the  al¬ 
gorithm.  If  each  point  keeps  q  records,  the  Bellman-Ford 
architecture  will  need  up  to  q{n  —  1)  rounds  of  relaxation 
operations,  so  a  small  q  is  better.  (In  our  implementation, 
we  choose  q  =  4.) 

The  bound  of  the  second  derivative  of  the  path.  A, 
enters  in  the  following  way.  Due  to  the  ^sumption  of 
constant  speed,  the  derivative  of  the  speed  comes  from  a 
direction  change.  Obviously,  a  piecewise  linear  trajectory 
has  unbounded  second-order  derivative.  We  therefore  limit 
the  impact  of  the  action  of  changing  directions:  the  action . 
doesn’t  finish  all  at  once;  instead,  it  takes  a  small  amount 
of  time  proportional  to  the  time  needed  to  traverse  distance 
A.  The  assumption  guarantees  that  there  is  enough  room 
to  finish  the  action  of  changing  directions.  Also,  fine 
discretization  reduces  the  capability  to  changing  directions 
from  one  segment  to  another.  As  \\p  =  1||,  the  angle  (f> 
between  two  consecutive  segments  is  limited  by 


Returning  to  step  2  of  the  relaxation  operation,  only  if  the 
angle  between  the  new  segment  i^k  and  the  last  segment  of 
the  path  pk  satisfies  this  inequality,  is  the  extension  feasible. 
The  corresponding  shift  is  called  a  feasible  shift. 

The  key  step  in  the  SSSP  algorithm  is  the  computation 
of  the  approximate  tilde  value  function  V{zk,Pk)  at  the 
shifted  point  Zk  in  step  1.  In  the  isotropic  case  this  can  be 
done  by  interpolation,  as  in  the  Fast  Marching  algorithm. 


Xjt+l 


Fig.  5.  Extension  of  a  shifted  path.  When  considering  the  next  segment 
TTfc  =  Pk  shifts  Zk  'rm  shape-preserving  way. 


What  we  want  here,  however,  is  a  path-oriented  method 
that  takes  as  an  input  only  the  shifted  path  Pk,Sk*  given 
by  the  recorded  discrete  path  pk  and  shifts  5^.  The  idea, 
roughly  speaking,  is  to  shift  the  path  in  a  “shape-preserving” 
way,  as  shown  in  Fig.  5.  When  the  end  point  of  the  path 
Pk  shifts  from  Xk  to  a  non-grid  point  Zk  =  Xk  -h 
the  trajectory  of  tiie  path  should  shrink  (or  swell)  in  such 
a  way  that  the  resulting  path  keeps  its  original  shape.  By 
doing  so,  the  directions  of  all  the  segments  of  the  path  are 
almost  the  same  as  those  before  the  shift.  We  call  such 
shifts  shape-preserving.  The  shape-preserving  shift  results 
in  a  backward  propagation  that  affects  the  shifts  of  all 
the  segments,  turning  shifts  Sk  =  (si,...,sjfc_i,0)  into 
Sk+i  =  (^1  +  Sij„.jSk^i  +  The  value  function 

over  the  shape-preserving  shifted  path  Pk,Sk+x^  is  thus  a  good 
approximation  of  the  tilde  value  function  V{zk,Pk)-  The 
collection  of  these  values  then  forms  a  profile  that  is  used 
to  compute  the  minimizing  shift  in  step  2. 

It  is  not  necessary  to  complete  the  backward  propaga¬ 
tion  before  computing  the  value  functions  of  the  shape¬ 
preserving  shifted  paths.  We  can  do  it  recursively  based  on 
the  knowledge  of  the  profile  of  V{zk^i,Pk-i)^  Let  sj^_j 
be  the  first-step  propagation  due  to  the  shape-preserving 
shift  sj..  We  have  Zk^i  =  Xfc_i  H-  (sjt-i  -f  and 

Zk-Xk+  s^kCk^  Thus, 

V{Zk,Pk)  =  V{Zk-l,Pk^lH\\^h-Zk^l\lc{XkyZk-Zk,^l)^ 

(21) 

Based  on  the  property  of  the  shape-preserving  shift,  the 
propagated  shift  depends  on  in  one  of  two  ways. 
The  first  is  the  proportional  shift.  It  occurs  when  segments 
TTk-i  and  TTk  have  the  same  baseline  direction  and  the 
same  quadrantal  direction.  With  proportional  shift,  = 
It  is  the  consequence  of  the  fact  that  x^,  the 
start  point  of  the  first  segment,  never  shifts.  If  the  path  turns, 
either  by  changing  the  baseline  direction  or  the  quadrantal 
direction,  the  propagated  shift  happens  in  a  parallel  manner. 
In  this  case  the  last  segment  tt^^i  retains  its  direction  and 
is  computed  accordingly. 

With  the  profile  of  the  approximate  tilde  value  function 
V{zkiPk)y  the  optimization  problem  in  step  2  becomes  a 
minimizing  problem  over  a  single  scalar  variable  s^,  thus 
is  significantly  simplified.  The  profile  is  in  a  piecewise 
linear  form  over  the  inteval  sjj.  6  H  =  [— 1»1]'  However, 
not  every  e  [—1,1]  is  feasible  even  if  it  satisfies  the 
direction  condition  (20).  In  the  example  shown  in  Fig.  6, 
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Fig.  6.  /B^  bounds  of  the  shifts.  The  dotted  lines  show  how  the 

path  extends  in  the  gird  ^ph;  the  dashed  lines  show  the  shift  paths;  and 
the  brackets  show  the  Bf/B]~  bounds  of  the  shifts,  ta  the  upper  figure, 
spans  the  (-1, 1]  interval.  When  the  path  extends  to  xs,  as  shovm 
in  the  lower  figure,  «2  bounded  by  B^  <1  from  above.  This  is  because 
the  shift  Si  of  the  extended  path  increases  in  the  shape-preserving  way, 
and  at  the  same  time  has  to  be  less  than  or  equal  to  1. 


Fig.  7.  This  figure  shows  how  SSSP  generates  a  loop  when  the  path 
extends  from  xi  to  X2,  then  back  to  X3  =  xi.  The  dotted  lines  show  how 
the  path  extends  in  the  grid  graph;  the  dashed  lines  are  the  shifted  paths. 
The  reason  of  generating  loops  is  that  the  tilde  value  function  decreases 
as  the  paths  extend,  as  shown  wi  the  right 


$2  is  bounded  by  9^  <1  because  si  after  the  backward 
propagation  must  be  within  its  [-1, 1]  interval.  We  denote 
these  bounds  9^/9^,  If  these  bounds  are  tight,  it  can  be 
proven  that  the  resulting  shifted  path  has  all  shifts  less  than 
or  equal  to  1. 

The  SSSP  algorithm  is  not  loop-free.  Although  the  cost 
function  c{x,  v)  is  strictly  positive  at  all  x  and  v,  the  tilde 
value  function  may  decrease  as  the  path  extends.  See  Fig.  7 
for  reference.  Suppose  the  (isotropic)  cost  function  takes 
only  two  values:  c(rr)  =  10  for  all  points  to  the  left  of  xi, 
and  c{x)  =  1  for  dl  points  to  the  right.  We  can  compute 
the  tilde  value  functions  V{xi,pi)  =  \/5  x  10  =  22.361, 
and  V^(a:2,P2)  =  2  x  10  4-  V5  =  21.414.  Therefore,  by 
extending  the  path  pi  to  rc2,  the  tilde  value  decreases  by 
0.947.  The  decreasing  tilde  value  implies  the  possibility  of 
loops.  In  fact,  remember  that  SSSP  records  multiple  paths  at 
each  point.  The  multiple  records  also  imply  potential  loops 
as  paths  extend.  We  definitely  don’t  want  looped  paths. 
One  way  of  preveting  loops  is  to  explicitly  check  the  paths 
and  delete  those  with  loops.  However,  the  looped  path  is 
inherent  within  the  shifted  path  context  and  makes  sense. 
In  the  lowest  part  of  Fig.  7,  when  we  extend  the  path  back 
to  X3  =  xi,  the  tilde  value  function  of  that  point  improves 
to  F(x3,P3)  =  2  X  10+>/5  =  22.236,  based  on  the  path  pz 
with  the  dashed  line.  Therefore,  though  ps  contains  loop, 
the  shifted  path  pz  doesn’t.  Fortunately,  pz  is  still  achievable 
by  replacing  xi  with  x^.  Thus,  the  explicit  loop  elimination 
doesn’t  affect  the  optimality  of  the  results. 

Unlike  random  sampling  algorithms,  which  use  up  run¬ 
ning  time  in  constructing  the  graph,  the  SSSP  algorithm 
spends  more  time  on  its  path  computation.  Specifically,  for 
a  graph  of  n  sample  points  and  m  links,  random  sampling 
algorithms  take  O(n^)  time  to  select  links  by  checking 


Fig.  8.  Path  conqjutation  problem  over  a  rectangular  region  with  flat 
cost  function  c(x,v)  =  0.1.  The  grid  in  dotted  line  shows  the  grid 
discretization  of  the  region.  Sfc*s  are  the  shifts.  The  triangle-labeled  discrete 
path  shows  how  serious  the  limit-on-direction  problem  could  be. 


the  distance  between  the  sample  points,  and  0{m)  time  to 
compute  the  values  of  each  link.  The  SSSP  algorithm,  on 
the  other  hand,  needs  0{qmn)  relaxation  operations,  where 
q  is  the  number  of  records  kept  at  each  point  of  the  graph. 
During  each  of  the  relaxation  operation,  SSSP  performs 
a  fixed  number  of  computations  to  form  a  profile  of  the 
tilde  value  function  K(z,p),  and  to  find  the  minimizing 
shift  s'.  The  time  complexity  of  the  SSSP  algorithm  is  thus 
0{qmn),  This  is  a  similar  expression  as  the  time  complexity 
of  the  random  sampling  algorithms.  However,  SSSP  can 
find  good  paths  even  with  a  coarse  grid  and  fewer  points 
in  the  graph.  Therefore,  the  running  time  of  SSSP  should 
be  smaller  than  that  of  the  random  sampling  algorithms,  as 
we  show  in  the  next  section. 
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Rg.  9.  The  typical  senario  and  minimum-risk  paths  computed  by  two 
algorithms,  (i).  the  proposed  SSSP  algorithm:  361  (19x19)  points,  path 
value  is  358.2,  running  time  is  2.6  sec;  (ii).  the  uniform  random  sampling 
algorithm:  3{X)0  points,  path  value  is  375.3,  running  time  is  74.6  sec. 


V.  Numerical  Results 

The  first  set  of  experiments  illustrates  the  effect  of 
the  shift  concept.  Fig.  8  shows  a  senario  where  the  cost 
function  is  0.1  at  all  positions  and  directions.  The  flat 
cost  function  represents  the  simplest  case  of  the  path 
computation  problem.  The  simple  example  demonstrates 
how  the  shiftment  helps  to  find  better  paths.  In  Fig.  8, 
an  8“by-8  region  is  discretized  with  grid  size  h  —  1,  We 
want  to  compute  the  shortest  path  from  the  starting  point 
Xs  =  (0,0)  to  the  destination  point  Xd  =  (8,6).  If  we 
run  Dijkstra’s  shortest  path  algorithm  over  the  grid  graph, 
the  triangle-labeled  discrete  path  is  the  solution  with  path 
value  being  1.4.  (Actually,  all  the  zigzag  paths  from  Xg 
to  Xd  are  equally  the  shortest  path.)  This  obviously  bad 
solution  can  be  improved  by  increasing  the  neighborhood 
distance,  which  allows  longer  segments  than  those  only 
between  the  adjacent  points  on  grids.  The  improved  discrete 
path  is  represented  by  the  dashed  line  and  diamond  label 
in  Fig.  8.  The  path  value  reduces  to  1.013.  The  discrete 
path  cannot  achieve  the  optimum  value  1.0  of  the  straight 
line  connecting  the  starting  and  destination  points,  because 
of  the  “limit-on-direction”  problem  mentioned  before.  By 
introducing  the  shifted  path,  the  segments  are  free  to  take 
different  directions,  thus  generate  the  near-to-optimal  path. 
The  shifted  path,  labeled  by  circles  with  path  value  1.0,  is 
very  close  to  the  optimal  straight  line.  We  also  display  in 
Fig.  8  the  amount  of  the  shifts  at  the  consecutive  grid  points 
of  the  discrete  path. 

We  now  apply  the  proposed  SSSP  algorithm  to  the  com¬ 
putation  of  paths  for  Unmanned  Air  Vehicles  (UAVs)  that 
minimize  the  risk  of  detection  (and  destruction)  by  ground 
radars  of  Surface-to- Air  Missile  (SAM)  sites.  Following  [1], 
finding  the  minimum  risk  path  is  abstracted  as  the  ASP 
problem  with  cost  function 
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Fig.  10.  Results  value  functions  of  the  SSSP  algorithm  and  random 
sampling  algorithm  at  different  san^ling  rates. 


where  Zi  is  the  position  of  the  ith  radar,  and 
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is  the  risk  density  function,  representing  the  probability 
of  the  UAV  being  detected  by  the  ith  SAM  site.  This 
probability  depends  on  the  aircraft’s  Radar  Cross  Section 
(RCS),  which  is  a  measure  of  its  ability  to  reflect  radar 
signals  in  the  direction  of  the  radar  receiver.  In  (23),  the 
RCS  is  an  ellipse  with  semi-minor  axis  a  and  semi-major 
axis  6.  Unlike  most  civilian  aircrafts  whose  RCS  is  close 
to  a  circle,  the  UAVs  have  flat  RCS  with  large  r  =  6/a. 
Therefore,  if  the  angle  between  the  moving  direction  v 
and  the  SAM  site  direction  (x  -  Zi),  7}{x,  v,  Zi)  takes  the 
small  value  a  if  ^  =  0  and  the  large  value  6  when  ^  =  90 
degrees. 

Fig.  9  considers  the  senario  with  two  long  range  SAM 
sites  and  one  short  range  SAM.  The  circles  show  the 
effective  ranges  of  the  SAM  radars.  The  figure  also  shows 
the  minimum-risk  paths  obtained  by  the  SSSP  algorithm  and 
the  uniform  random  sampling  algorithm,  respectively.  The 
SSSP  algorithm  runs  over  a  19x  19  grid,  or  361  points  and 
yields  a  shifted  path  with  tilde  value  358.2.  The  uniform 
random  sampling  algorithm,  at  the  same  time,  runs  over 
3000  randomly  sampled  points  and  yields  a  discrete  path  of 
value  375,3.  Moreover,  the  running  time  of  latter  is  nearly 
30  times  slower  than  that  of  SSSP.  The  next  two  figures. 
Fig.  10  and  11,  show  more  results  of  the  two  algorithms  as 
we  change  the  sampling  rates.  The  SSSP  algorithm  keeps 
outperforming  the  uniform  random  sampling  algorithm  on 
both  the  path  values  and  the  running  time.  They  also  demon¬ 
strate  that  SSSP  is  able  to  achieve  the  optimal  solution  to 
the  ASP  problem  with  a  very  low  sampling  rate. 


VI.  Conclusion 

We  studied  the  problem  of  computing  the  shortest  paths 
over  a  continuous  region  in  this  paper.  The  cost  depends 
on  both  the  position  and  direction  of  movement,  thus  form¬ 
ing  an  Anisotropic  Shortest  Pathing  (ASP)  problem.  We 
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Fig.  11.  Results  running  time  of  the  SSSP  algorithm  and  random  sampling 
algorithm  at  different  sampling  rates. 

proposed  a  Shiftable  Segment-by-Segment  Pathing  (SSSP) 
algorithm  that  has  the  following  features. 

First,  SSSP  uses  rectangular  grid  to  discretize  the  region. 
Compared  with  other  sampling  schemes,  the  grid  graph  is 
easy  to  construct  and  is  independent  of  the  cost  function. 
Second,  the  SSSP  algorithm  uses  a  Bellman-F6rd  architec¬ 
ture  to  perform  a  breadth-first  search.  Each  point  of  the 
graph  keeps  multiple  records  of  most  promising  paths  to 
assure  the  optimality  of  the  result.  The  third  feature  is  the 
introduction  of  shifted  paths,  which  allow  the  discrete  paths 
to  shift  away  from  the  grid  points,  thereby  overcoming  the 
limit-on-direction  problem  that  arises  in  grid  discretization. 
The  last  feature  is  the  heuristic  to  compute  the  tilde  value 
function  V{x^p\,  By  constructing  the  profile  of  the  tilde 
value  functions  V (z,p)  over  the  shifted  points  z,  we  replace 
the  optimization  problem  over  all  shifts  to  the  simpler  one 
over  only  the  last  shift  of  the  path.  The  key  technique  of  the 
heuristic  is  the  shape-preserving  shift,  which  decides  how 
the  shifted  path  changes  when  the  end  point  shifts. 

These  features  of  SSSP  woric  together  and  lead  to  an 
efficient  algorithm.  The  theoretical  analysis  says  that  the 
SSSP  algorithm  has  a  similar  time  complexity  as  the  ran¬ 
dom  sampling  algorithms.  The  numerical  experiments,  on 
the  other  hand,  show  that  SSSP  outperforms  the  uniform 
random  sampling  algorithm  on  both  the  path  values  and  the 
running  time. 

Although  we  described  and  implemented  the  SSSP  algo¬ 
rithm  for  two-dimensional  regions,  it  is  easy  to  generalize 
the  algorithm  to  higher  dimensional  spaces.  The  major 
modifications  focus  on  the  local  shift  area  (I,  the  shifts  SkS 
and  the  computation  of  the  profile  of  V{zyp),  For  example, 
in  a  three-dimensional  space,  the  shift  is  restricted  within  a 
small  two  dimensional  plane,  ft(ar),  around  the  grid  point 
X,  Of  course,  the  algoritfim  will  ne^  more  computations  to 
construct  the  profile  of  V{zyp)  over  n(x). 
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